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An algorithm is presented which generates pairs of oscillatory random time 
series which have identical periodograms but differ in the number of oscillations. 
This result indicate the intrinsic limitations of spectral methods when it comes 
to the task of measuring frequencies. Other examples, one from medicine and 
one from bifurcation theory, are given, which also exhibit these limitations of 
spectral methods. For two methods of spectral estimation it is verified that the 
particular way end points are treated, which is specific to each method, is, for 
long enough time series, not relevant for the main result. 
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1 Introduction 



Nearly everywhere in science we encounter the problem of determining the fre- 
quency of oscillation of some autonomously oscillating system from which a 
signal has been measured over a finite time period. Most commonly, spectral 
methods are \iscd: From the squared modulus of the Fourier transform of the sig- 
nal (the periodogram) its expectation value (the power spectrum) is estimated, 
and the frequency is determined by some characterization of the position of a 
peak in the power spectrum (see, e.g., Mo et ai, 1988; Muzi & Ebert, 1993; 
Timmer et ai, 1996; Mcndez et ai, 1998; Godano & Capuano, 1999; Korhonen 
et ai, 2002; Slavic et ai, 2002). For linear systems driven by white noise it is 
easily seen that the periodogram or, equivalently, the empirical autocorrelation 
functicm do, in fact, contain all relevant information (Brockwell & Davis, 1991; 
Priestley 1981). 

For nonlinear oscillators, this is generally not the case. As was shown re- 
cently (Rossberg, 2002), the power spectrum of a noisy or chaotic (i.e., not 
perfectly periodic) nonlinear oscillator which oscillates with a given average fre- 
quency , i.e. the number of cycles per unit time, can be completely arbitrary. 
In particular, the power spectrum is unrelated to the frequency. The number of 
oscillation cycles is here defined by the number of transition of the trajectory of 
the signal in delay space through a Poincare section. Since this corresponds to 
a purely topological relation between the trajectory and the Poincare section, 
the corresponding measure of frequency was named topological frequency cjtop- 
Of course, intuition and experience suggest that, for typical experimental time 
series, the result of using some standard method to determine the oscillation 
frequency from the power spectrum would not differ too much from the topo- 
logical frequency. But from the power spectrum alone no upper limit for this 
difference can be derived. Nor seems any other efficient method to be known 
which would allow to determine an upper limit for this difference. 

After giving some examples for oscillations which exhibit this unrelatedness 
of topological frequency and power spectrum in Section 2, a result is presented 
which is in some sense stronger than the result cited above. A precise determi- 
nation of frequency is not only impossible from the power spectrum of a time 
series, but also from the periodogram. Since the power spectrum is, in one 
way or another, estimated from the periodogram by an averaging procedure, 
the estimated power spectrum contains less information than the periodogram. 
One might conjecture that from the additional information contained in the 
periodogram, such as higher order correlations between the amplitudes of in- 
dividual Fourier modes, frequency information can be extracted. As is shown 
in Section 3, this is not the case: Examples for time series are given which 
obviously have different topological frequencies, but, by construction, identi- 
cal periodograms. The use of the periodograms of the raw data is justified for 
spectral estimation only if end effects can be neglected, for example, when the 
time series are long enough. In Section 4 it is shown that, in fact, for long 
enough time series the result holds valid also when end effects are explicitly 
taken account by spectral estimators. 
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2 Oscillation frequencies and power spectra of 
some nonlinear oscillators 



As a first example, consider the noisy, weakly-nonlinear oscillator described by 
a complex amplitude A{t) with dynamics given by the noisy Landau-Stuart 
equation (Risken, 1989) 

i = (e + iuJo)A - (1 + ^g^)\A\^A + r){t), (1) 

where e, uq, and gi are real and r]{t) denotes complex, white noise with correla- 
tions 

{rj{t)rj{t'))^0, {rj{tMtr)^m-t') (2) 

[* = complex conjugation, (•) = expectation value]. In a certain sense (Arnold, 
1998), this system universally describes noisy oscillations in the vicinity of a 
Hopf bifurcation. In general (i.e., with gi =/= 0) the linear frequency loq, the 
spectral peak frequency Wpeak, i-e. the position of the maximum of the power 
spectrum Sa{co), the average frequency or phase frequency 

Wph,A := (tOi) , := Im{y4/A} , (3) 

and the mean frequency of A{t) 

_ (u^ilAI^) _ Im {AA*) _ JujSAito)dLU 

are all different; see Fig. 1. [Definitions (3,4) are sometimes restricted to an- 
alytic signals, characterized by 5a(w) = for w < 0, which are in one-to-one 
correspondence with the zero-mean, real- valued signals Re{A{t)}. See Boashash 
(1992) for the history.] 

The phase frequency measures the average number of circulations around 
the point A = in phase space per imit time (decompose A{t) = a{t)e^'^^*'> 
to see this). Thus, with the natural assumption that a circulation around the 
origin of phase space corresponds to a cycle of oscillation - which is justified by 
symmetry considerations - the phase frequency measures cycles per unit time, 
i.e., the oscillation frequency. Although this concept of an oscillation cycle 
is slightly different from that used for the definition of topological frequency, 
both frequency measures often give identical results, in particular when the 
oscillation amplitude does not come too close to zero. For a discussion of the 
effect of filtering on the phase frequency, see Rossberg (2002). In contrast to the 
phase frequency, the spectral peak frequency uJpcak.A and the mean frequency 
Wmean,^ are frequency measures which can be obtained directly from the power 
spectrum. As can be seen from Eqs. (3) and (4), Wph,A = Wmean.A if and only 
if and \A\'^ are uncorrelated. This is always the case for linear, noise-driven 
oscillators, but for general nonlinear oscillators it is not. 

When it is known that the dynamics of A{t) is of the form of Eqs. (1) and 
(2) and the power spectrum Sa{<^) of A{t) is precisely determined, then it is of 
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course possible to obtain the free parameters of Eq. (1) from a fitting procedure 
and to use these to determine u)ph,A- But none of these assumptions is typically 
satisfied in practice. Equations (1) and (2) are exact only in a particular limit, 
and 5/1 (w) is usually biased by some nontrivial transfer function inherent in the 
measurement process. Then the information about Wph,A cannot be extracted 
from the power spectrum anymore. The situation becomes aggravated by the 
fact that, as was shown by Scybold & Risken (1974), S'a(w) is always well 
approximated by a Lorenzian, uniformly for all e. 

As a second example, an experimental time series from the hand tremor of 
a subject diagnosed with Parkinson's disease shall be examined. The data was 
recorded by an accelerometer attached to the hand. It consists of 10240 samples 
recorded at a rate of 300 Hz. The data is available on the internet^. 

Figure 2a shows part of the time scries. In order to obtain an estimate of 
the power spectrum (Fig. 2b), the periodogram of the (tapered) time series was 
smoothed using a window of width h= 10 with weights w{i) = max(0, l—i^/h?). 
Its peak frequency is at (jjpcak/27r = 9.49 Hz, corresponding to 324 oscillations 
over the time series. In order to determine the phase frequency, the (band- 
pass filtered) analytic signal corresponding to the time series is obtained by 
convoluting it with the Morlet wavelet 



with loq = cijpoak and Ao; = cijpoak/6. The resulting phase frequency Uph = 
9.58 Hz (Fig. 2b) is robust under large variations of the filter width Aw and 
mid frequency ujq. It is also identical to the topological frequency obtained 
by band-pass filtering the time scries by a convolution with Re{w(t)} and a 
successive 2D delay embedding, e.g., with delay r = (23/300) s. Unfortunately, 
the investigated time series is too short to yield the observed difference between 
peak frequency and topological frequency significant in a statistical sense, as a 
careful statistical analysis using the method of Timmer et al. (1997) shows: The 
95% confidence interval for the peak frequency extends from 9.40 Hz to 9.61 Hz. 
It would be desirable to repeat the analysis on longer samples. 

It should be emphasized that the point of Rossberg (2002) is that the topo- 
logical or phase frequency cannot be derived from the power spectrum as a 
whole. Determination of the spectral peak frequency is only the most obvious 
attempt to do so. The failure of this particular approach has been noticed 
before, see, e.g., Pikovsky et al. (1997). 

3 Does the periodogram help? 

Is it possible to identify a relative shift between spectral peak frequency and 

topological or phase frequency by using a different method to evaluate the peri- 
odogram? As mentioned before, the answer to be given here is negative: there is 

'^http : //webber . physik . iiiii-f reiburg . de/^ j et i/tremordaten/park/0308-2-ali . dat 
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no general way to tell from the periodograms (or the empirical autocorrelation 
functions) of time series that they oscillate with different frequencies. 

For simplicity, only real valued time series x{t) shall be considered here for 
which the number of oscillations in a given time interval of length T can be de- 
termined by a simple 2D delay embedding. If there is a delay At, such that the 
trajectory of the delay vector {x{t),x{t — At)) encircles the origin of the coor- 
dinate system at finite distance, the number of circulations m gives the number 
of oscillations in the time series. The resulting frequency is then 2'Km/T . This 
is the definition of topological frequency restricted to 2D embedding. It is very 
similar to the one that can be obtained from of the definition of instantaneous 
phase (t){t) by Pikovsky et al. (1997) when defining u as the temporal average 
of d(j)/dt. However, their definition is based on the actual phase space of a dy- 
namical system and not on its approximate reconstruction from a time series by 
delay embedding. This distinction becomes particularly relevant when random 
time series, as they are modeled by the method described here, are evaluated. 

We shall now outline an iterative algorithm that randomly generates pairs of 
time series Xi{t), X2{t) of the type described above. The two time series gener- 
ated by this algorithm have identical periodograms but differ in the number of 
oscillation cycles. (The algorithm can also be considered as a source of random 
time series Xi(t), X2{t) which have identical power spectra but different expec- 
tation values for the oscillation frequencies - in fact, the frequencies are fixed). 
The algorithm uses some ideas from an algorithm of Schreiber & Schmitz (1996) 
that generates time series with periodograms nearly identical and distributions 
of values identical to a given times series. 

The main part of our algorithm generates two complex- valued random time 
series yi{t), 2/2 (t) with N samples at t = 0, . . . , iV — 1. When continued peri- 
odically [yi{t) = yi{t + N)], the distribution of the output is invariant under 
translation in time. The periodograms of the two time series are identical. The 
moduli \yi{t)\ follow a predetermined narrow distribution with values near one, 
i.e., the values of yi{t) lie near the unit circle in the complex plane. yi{t) and 
y2{t) encircle the origin of the complex plane ki and k2 (fci ^ ^2) times re- 
spectively. Using these time series, the oscillatory, real- valued time series xi{t), 
X2{t) {xi{t) = 'Re{exp{iuiot)yi{t)) with identical periodograms are constructed. 

In order to obtain the final time series yi (t) , the following two basic steps are 
performed iterative, starting with two time series which encircle the origin of the 
complex plane fci respectively k2 {ki ^ ^2) times, but do not yet have identical 
periodograms: In Fourier space, the magnitudes of the Fourier modes of the 
two time series are adjusted, such that the magnitudes of the Fourier modes 
at the same frequencies approach each other, while the phases of the Fourier 
modes are kept constant (Steps 4 and 5 in the detailed description below). 
Then, back in physical "space", i.e., using the actual time series, a nonlinear, 
monotonously increasing transformation of the magnitudes of the values for each 
time is applied in such a way that the desired distribution of the complex time 
series near the unit circle (see above) is enforced (Step 10 below). If the phases 
were kept constant in this second basic step, the algorithm would (as numerical 
experiments show) typically converge to a trivial solution with yi{t) = 2/2 (t)- 
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Obviously, the number of circulations of the origin of the complex plane would 
then have changed with respect to the initial condition - at least for one of 
the two time series. In order to avoid these phase slips, the difference in phase 
angle between successive values is, for both time series, enforced to lie within 
a certain, small rage, by adjusting, at each iteration, also the phase angles 
(Steps 1,6-9). The details of these manipulations of the phases, as well as the 
precise choice of the initial time scries (Steps 2 and 3), contain random elements, 
as will be explained below. Experience shows that convergence of the algorithm 
is, although not guaranteed, with appropriate choices of the tuning parameters 
of the algorithm (see below), highly probable. 

This is a step-by-step description of this algorithm. We give the values of 
tuning parameters as they were used to generate the sample Pair A shown in 
Fig. (3), which consists of two time-series with N = 256 points, in parenthesis. 

1. Generate two random permutations : (0 . . . (A^ — 1)) ^ (0 . . . (A^ — 1)) 
(/ = 1, 2) of length N, with equal probability 1/iV! for each permutation, 
to be used later. 

2. Initialize the series yi{t) and 2/2 (i) as yi{t) = exTp{27rikit/N) a.t t = 
0, . . . ,N — 1 with small integer fc; (e.g., fci = and k2 = 1). 

3. In order to excite also the remaining Fourier modes, add some noise to 
yi{t) and y2{t) (c-g-, complex Gaussian with variance 2 x 10~^). 

Now the main loop (steps 4-10) starts. 

4. Calculate the discrete Fourier transforms yiicu) of the time series yi{t). 
Increase or decrease the values of \yi{to)\ such that the distance between 
|yi(a;)| and |if2('i^)| decreases by a fixed percentage (e.g., 20%) while keep- 
ing arg(y;(a;)) and \yi{Lo)\ + \y2{Lo)\ fixed. In this step the two periodograms 
are forced to approach each other. Notice, however, that the manipula- 
tions in steps 6-10 may have an adversary effect. In order to have the 
algorithm as a whole converge, some balancing of the tuning parameters 
is therefore required. 

5. Obtain the inverse Fourier transforms of the results and assign them to 
yi{t) and y2{t). 

6. Calculate the increments of a.rg{yi{t)) between successive samples. Apply 
the permutations 11; to the sequences of phase increments (see below for 
an explanation). 

7. For phase increments of absolute value larger than some threshold (e.g., 
0.4) distribute part of the increment (e.g., 1%) symmetrically over the pre- 
ceeding and succeeding increments such that the sum of the three remains 
unchanged. Do this step with parallel update for all increments. 

8. Repeat step 7 until no phase increment exceeds the threshold. 
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9. Obtain the inverse permutations of these phase increments, and sum 
them up to obtain two new time series for the phases aTg{yi{t)). 

10. For each series, find a polynomial of low order (e.g., 9th order), such 
that the distribution of Pi{\yi{t)\) fits best with the desired distribution 
of the arguments (we use a Gaussian with mean m = 1 and variance 

= 0.04).^'^ Construct two new complex-valued time series by using the 
arguments Pi{\yi{t)\) and the phases from step 9, and assign them to yi{t) 
and y2{t). 

11. Repeat steps 4-10 until the periodograms of yi{t) and yi{t) converge to 
the same values, i.e., 6 ^ "'^ X^wdyil'^)! ^ |2/2('jj)|)^ < e, with some 
small e > (e.g., e = 10^^). For larger N convergence is substantially 
improved by judiciously choosing new permutations 11; from time to time 
(see below). 

12. Multiply random phases exp(i(f>i) to yi{t) and y2{t) in order to make the 
distributions of xi{t) and X2{t) obtained in the next step time-translation 
invariant. 

13. Calculate the real- valued time series xi{t) — Rc{exp{iu!ot)yi{t)} with some 
not too small wq (e.g., ivq = 57r/32). 

14. As in steps 4 and 5, remove the small differences in the periodograms of 
Xi{t) and X2(t) resulting from the Re{-} operation in step 13. But now, 
do it in a single (100%) step. 

A typical result of running this algorithm is shown in Fig. 3. With a high enough 
frequency of the fast modulation u)q, the algorithm guarantees that the origin 
of delay space [At « 7r/(2cjo)] is encircled at finite distance (see Fig. 4). 

Now, some details concerning the manipulation of phases (Steps 6-9) shall 
be explained. The permutations of the phase increments prior to the relaxation 
of large phase jumps in steps 7 and 8 prevent the algorithm from running into 
solutions with localized regions of fast drift of phase, which do not occur in 
typical experimental time series. However, it turns out that, with permutations 

^Pj(-) is determined such that 

JV-l 
min= J2 
s=0 

where Qi{-) is the permutation which rank orders the sequence {|j/i(t)|}t and erf~^(-) denotes 
the inverse error function given by 

X = —= / e " du for all -1 < a; < 1. (7) 

■'Schrcibcr & Schmitz (1996) achieve a similar result by using, instead of P;(-)i functions 
M(-) which map the sets of values of the time scries (c.g.,{|3/i(t)|}t) onto a given target set of 
values in such a way that rank ordering is preserved. The method used here has the advantage 
over their method that the periodograms can be matched to arbitrary accuracy. 



Pi {\yi (QiW)l) ■ 



■2l/2aerf- 



2s + 1 
N 



(6) 
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selected with equal probability from the set of all permutations, the algorithm 
has difficulties to converge for large N, while without any permutation it readily 
does. This dilemma can be solved by slowly interpolating from fully random 
permutations to the identity permutation: Execution of the main loop (steps 
4-10) is dived into blocks of, e.g., 200 iterations. After terminating a block, 
the yi{t) are reset to the pair of series with the lowest 6 so far, and the ran- 
dom permutations Hi are replace by the permutations which rank order the 
sequences {{t + hvi{t)) mod A^}f (t = Q, . . . , N — 1). Here {i'i{t)}t are sequences 
of independent random numbers distributed equally between and 1 which are 
renewed after each block and, as usual, uraoAN := u — mN, where m is the 
largest integer with u — mN > 0. The parameter h controls the interpolation 
between fully random permutations (h — > oo) and the identity (h = 0). Succes- 
sively smaller values of h are used after each block. As an example for a pair 
of long time series (iV = 8192, —ki = = 24), Pair B was generated using 
h = 400, 200, 100, 0, (convergence is reached after 5 blocks) and all other pa- 
rameters as for Pair A. For these long time series the fluctuations of the phases 
are clearly superimposed by a linear drift of the phase difference (Fig. 5). For 
Pair A this is not obvious, and one might be tempted to attribute the difference 
in frequency there to a "mere fluctuation" . A Matlab code that implements the 
algorithm in the form used to generate Pair B is available on the internet^. 

4 The influence of end-point effects for two spec- 
tral estimators 

Calculating the periodogram of a time series (and the empirical autocorrelation 
function as the Fourier transform thereof) is the correct first step toward esti- 
mating its power spectrum only if end-point effects can be neglected. This is 
the case when the time series allows a meaningful periodic continuation (which 
the algorithm presented here guarantees), or if the time series are sufficiently 
long. For two methods of spectral estimation it shall here be demonstrated 
that the condition of sufficient length can also be met by the generated pairs of 
time series - without making used of periodic continuation. That is, although 
the two methods differ in the way endpoints are treated, in the way the pe- 
riodograms are averaged, and also slightly in their results, each of them gives 
virtually identical results for the two time series of a pair. 

First consider Welch's method (Welch, 1967). Both time series of Pair B are 
split into 126 segments of 192 samples length with 2/3 overlap. The segments 
are tapered using a Hamming window and, in order to interpolate the resulting 
power spectra, extended to length 2048 by appending zeros. The periodograms 
of all segments of a time series are averaged. As is shown in Fig. 6, the estimated 
power spectra are nearly identical for both series. The differences are much 
smaller than then the difference of the oscillation frequencies. 

Next, autoregressive spectral estimation is used. Applying Burg's method 

^http : //www . f dm . uni-f reiburg . de/~axel/onthelimit s/ 
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(Marple, 1987), a p-th order autoregressive process 



x{t) = akx{t - k)j + e{t) (8) 

with uncorrelated, identically distributed e{t), is fitted to the two time series 
of Pair B. The resulting estimates for the coefficients (ai, .... dp) for p = 10 are 
listed in Table 1. No attempt was made to find the "optimal" model order p, 
since for all reasonable values of p the same observation is made: The power 
spectra of the two fitted autoregressive processes are virtually indistinguishable. 
For example, with p = 10 the spectral maxima arc both located at the angular 
frequency oj := 27r/ = 0.4887. They differ only by Au = 2 x 10~^. In order to 
get an estimate for the errors of the autoregressive spectral estimates, the first 
and second half of the two time series are fitted separately. Then the enforced 
property of having equal periodograms is lost and the estimated spectra differ 
slightly. The maxima of the resulting power spectra are located at w = 0.4876 
and 0.4897 for xi{t), and 0.4890 and 0.4884 for X2{t). For comparison, the 
oscillation frequencies of xi{t) and X2{t) are co = 0.4725 and 0.5093 respectively 
(see inset of Fig. 6). Hence, it is not possible to explain the undetectability 
of the difference in frequency between xi{t) and X2{t) of Pair S by a lack of 
spectral resolution. 



5 Conclusion 

Several examples were given that demonstrate that the power spectra of time 
series alone do not contain sufficient information to ensure that the oscillation 
frequencies (cycles per unit time) of the time series could precisely be deter- 
mined. These examples illustrate a corresponding general theorem (Rossberg, 
2002). Furthermore, it was shown that this result does not depend on the partic- 
ular method of spectral estimation, as long as it is functional in the periodogram 
or the empirical autocorrelation function, respectively. There is also no other 
method of evaluating the periodogram which would yield the oscillation fre- 
quency. This is shown by introducing a general iterative method to generate 
pairs of time series of arbitrary length which have identical periodograms but 
different frequencies of oscillation. 
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Table 1: The autoregressive coefficients in Eq. (8) with p = 10 estimated for 
Pair B (see text) using Burg's method. The variance of e{t) is estimated as 
var(e) = 0.0409 for both time series. 
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Figure Captions 
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Figure 1: The power spectrum S'A(t^) of A{t) given by Eqs. (1,2) with e = 2 and 
gi = 1, obtained from a numerical simulation, and cjpoak.A [sec Seybold & Risken 
(1974) for analytic results], compared to the mean frequency Wmean.A = u>o — 
5i(|A|4)/(|^|2) [(|^|2n^ ^ 2"Ar-id"Ar/de", where AT := 77^2 exp(eV4)(l + 
erf(e/2)), see Risken (1989)] and the phase frequency u)ph,A = wq — gi{^\A\^, 
defined by Eqs. (3,4), and the linear frequency uq. 
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Figure 2: Power spectrum and oscillation frequency (b) of a human hand- tremor 
time scries (a). In part (b). the vertical line corresponds to the phase frequency 
of the time series, the horizontal bar shows the 95% confidence interval for the 
spectral peak frequency. 
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Figure 3: (a) A pair of random time series xi (t), X2{t) {Pair A) with frequencies 
/ = 20/256 and / = 21/256 respectively. It was generated using the algorithm 
described in this report with the proposed parameters, (b) The periodogram, 
which is identical for both time series. 
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Figure 4: The trajectory of xi{t) in delay space with At = 3. By the fact that 
it is winding 20 times around the origin at a finite distance, the frequency is 
sharply defined. 
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Figure 5: The phase angles of the complex time series yi{t) (lower) and y2{t) 
(upper) that were used to generate xi{t) and X2{t) of Pair B (sec text). The 
two time scries differ by 48 cycles, and so do the time scries Xi{t) and X2{t). 
Notice that yi{t) x cxp{iuJot) and y2(i) x cxpl^iujot) can be approximated by the 
analytic signals corresponding to Xi{t) and X2{t). 
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Figure 6: Estimated power-spectral densities (PSD. in miits of power per radians 
per sample) for Pair B. The two estimates (dash-dotted) using autoregressive 
modeling (see text) are indistinguishable. When using Welch's Method (see 
text), the differences between the estimates corresponding to xi{t) (solid) and 
X2{t) (dotted) are hardly visible. The blowup of the peak region in the inset also 
shows the oscillation frequencies of xi(t) and X2{t) (vertical dashed lines) and, 
as an estimate of the nominal error of the autoregressive spectral estimates, the 
interval covered by the maxima of the spectra estimated by using only half the 
data of each time series (vertical solid lines near u) axis). Notice that the inset 
has all-linear axes, while the large graph is semi-logarithmic. 
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